home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
lin_alg.lha
/
lin_alg
/
vvector.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-08
|
11KB
|
378 lines
// This may look like C code, but it is really -*- C++ -*-
/*
************************************************************************
*
* Verify Operations on Vectors
*
************************************************************************
*/
#include "LinAlg.h"
#include <math.h>
#include <builtin.h>
#include <ostream.h>
/*
*------------------------------------------------------------------------
* Service validation functions
*/
static void verify_element_value(const Vector& v,const REAL val)
{
register imax = 0;
register double max_dev = 0;
register int i;
for(i=v.q_lwb(); i<=v.q_upb(); i++)
{
register double dev = abs(v(i)-val);
if( dev >= max_dev )
imax = i, max_dev = dev;
}
if( max_dev == 0 )
return;
else if( max_dev < 1e-5 )
message("Element #%d with value %g differs the most from what\n"
"was expected, %g, though the deviation %g is small\n",
imax,v(imax),val,max_dev);
else
_error("A significant difference from the expected value %g\n"
"encountered for element #%d with value %g",
val,imax,v(imax));
}
static void verify_vector_identity(const Vector& v1, const Vector& v2)
{
register imax = 0;
register double max_dev = 0;
register int i;
are_compatible(v1,v2);
for(i=v1.q_lwb(); i<=v1.q_upb(); i++)
{
register double dev = abs(v1(i)-v2(i));
if( dev >= max_dev )
imax = i, max_dev = dev;
}
if( max_dev == 0 )
return;
if( max_dev < 1e-5 )
message("Two #%d elements of the vectors with values %g and %g\n"
"differ the most, though the deviation %g is small\n",
imax,v1(imax),v2(imax),max_dev);
else
_error("A significant difference between the vectors encountered\n"
"at #%d element, with values %g and %g",
imax,v1(imax),v2(imax));
}
/*
*------------------------------------------------------------------------
* Test allocation functions and compatibility check
*/
static void test_allocation(void)
{
cout << "\n\n---> Test allocation and compatibility check\n";
Vector v1(20);
Vector v2(1,20);
Vector v3(0,19);
Vector v4(v1);
cout << "The following vector have been allocated\n";
v1.info(), v2.info(), v3.info(), v4.info();
cout << "\nStatus information reported for vector v3:\n";
cout << " Lower bound ... " << v3.q_lwb() << "\n";
cout << " Upper bound ... " << v3.q_upb() << "\n";
cout << " No. of elements " << v3.q_no_elems() << "\n";
cout << " Name " << v3.q_name() << "\n";
cout << "\nCheck vectors 1 & 2 for compatibility\n";
are_compatible(v1,v2);
cout << "Check vectors 1 & 4 for compatibility\n";
are_compatible(v1,v4);
#if 0
cout << "Check vectors 1 & 3 for compatibility\n";
are_compatible(v1,v3);
#endif
cout << "v2 has to be compatible with v3 after resizing to v3\n";
v2.resize_to(v3);
are_compatible(v2,v3);
Vector v5(v1.q_upb()+5); v5.set_name("Vector v5");
cout << "v1 has to be compatible with v5 after resizing to v5.upb\n";
v5.info();
v1.resize_to(v5.q_no_elems());
are_compatible(v1,v5);
cout << "\nDone\n";
}
/*
*------------------------------------------------------------------------
* Test uniform element operations
*/
static void test_element_op(const int vsize)
{
const double pattern = PI;
register int i;
cout << "\n\n---> Test operations that treat each element uniformly\n";
Vector v(-1,vsize-2);
Vector v1(v);
cout << "Writing zeros to v...\n";
for(i=v.q_lwb(); i<=v.q_upb(); i++)
v(i) = 0;
verify_element_value(v,0);
cout << "Clearing v1 ...\n";
v1.clear();
verify_element_value(v1,0);
cout << "Comparing v1 with 0 ...\n";
assure(v1 == 0, "v1 is not zero!");
cout << "Writing a pattern " << pattern << " by assigning to v(i)...\n";
for(i=v.q_lwb(); i<=v.q_upb(); i++)
v(i) = pattern;
verify_element_value(v,pattern);
cout << "Writing the pattern by assigning to v1 as a whole ...\n";
v1 = pattern;
verify_element_value(v1,pattern);
cout << "Comparing v and v1 ...\n";
assure(v == v1, "v and v1 are unexpectedly different!");
cout << "Comparing (v=0) and v1 ...\n";
assert(!(v.clear() == v1));
cout << "\nClear v and add the pattern\n";
v.clear();
v += pattern;
verify_element_value(v,pattern);
cout << " add the doubled pattern with the negative sign\n";
v += -2*pattern;
verify_element_value(v,-pattern);
cout << " subtract the trippled pattern with the negative sign\n";
v -= -3*pattern;
verify_element_value(v,2*pattern);
cout << "\nVerify comparison operations\n";
v = pattern;
cout << " (v=pattern)>0\n";
assure( v > 0, "Unexpectedly, v is not greater than 0");
cout << " (v=pattern)>-pattern\n";
assert( v > -pattern );
cout << " (v=-pattern)<-pattern/2\n";
v -= 2*pattern;
assert( v < -pattern/2 );
cout << "\nAssign 2*pattern to v by repeating additions\n";
v = 0; v += pattern; v += pattern;
cout << "Assign 2*pattern to v1 by multiplying by two \n";
v1 = pattern; v1 *= 2;
verify_element_value(v1,2*pattern);
assert( v == v1 );
cout << "Multiply v1 by one half returning it to the 1*pattern\n";
v1 *= 1/2.;
verify_element_value(v1,pattern);
cout << "\nAssign -pattern to v and v1\n";
v.clear(); v -= pattern; v1 = -pattern;
verify_element_value(v,-pattern);
assert( v == v1 );
cout << "v = sqrt(sqr(v)); v1 = abs(v1); Now v and v1 have to be the same\n";
v.sqr();
verify_element_value(v,pattern*pattern);
v.sqrt();
verify_element_value(v,pattern);
v1.abs();
verify_element_value(v1,pattern);
assert( v == v1 );
cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1\n";
for(i=v.q_lwb(); i<=v.q_upb(); i++)
v(i) = 2*PI/v.q_no_elems() * i;
v1 = v;
v.user_function(sin);
v1.user_function(cos);
v.sqr();
v1.sqr();
v += v1;
verify_element_value(v,1);
cout << "\nVerify constructor with initialization\n";
Vector vi(1,5,0.0,1.0,2.0,3.0,4.0,"END");
Vector vit(5);
for(i=vit.q_lwb(); i<=vit.q_upb(); i++)
vit(i) = i-1;
verify_vector_identity(vi,vit);
cout << "\nDone\n\n";
}
/*
*------------------------------------------------------------------------
* Test binary vector operations
*/
static void test_binary_op(const int vsize)
{
const double pattern = PI;
register int i;
cout << "\n---> Test Binary Vector operations\n";
Vector v(2,vsize+1);
Vector v1(v);
Vector vp(v);
for(i=v.q_lwb(); i<=v.q_upb(); i++)
vp(i) = (i-v.q_no_elems()/2.)*pattern;
cout << "\nVerify assignment of a vector to the vector\n";
v = pattern;
v1.clear();
v1 = v;
verify_element_value(v1,pattern);
assert( v1 == v );
cout << "\nAdding one vector to itself, uniform pattern " << pattern << "\n";
v.clear(); v = pattern;
v1 = v; v1 += v1;
verify_element_value(v1,2*pattern);
cout << " subtracting two vectors ...\n";
v1 -= v;
verify_element_value(v1,pattern);
cout << " subtracting the vector from itself\n";
v1 -= v1;
verify_element_value(v1,0);
cout << " adding two vectors together\n";
v1 += v;
verify_element_value(v1,pattern);
cout << "\nArithmetic operations on vectors with not the same elements\n";
cout << " adding vp to the zero vector...\n";
v.clear(); v += vp;
// verify_vector_identity(v,vp);
assert( v == vp );
v1 = v;
cout << " making v = 3*vp and v1 = 3*vp, via add() and succesive mult\n";
add(v,2,vp);
v1 += v1; v1 += vp;
verify_vector_identity(v,v1);
cout << " clear both v and v1, by subtracting from itself and via add()\n";
v1 -= v1;
add(v,-3,vp);
verify_vector_identity(v,v1);
cout << "\nTesting element-by-element multiplications and divisions\n";
cout << " squaring each element with sqr() and via multiplication\n";
v = vp; v1 = vp;
v.sqr();
element_mult(v1,v1);
verify_vector_identity(v,v1);
cout << " compare (v = pattern^2)/pattern with pattern\n";
v = pattern; v1 = pattern;
v.sqr();
element_div(v,v1);
verify_element_value(v,pattern);
compare(v1,v,"Original vector and vector after squaring and dividing");
cout << "\nDone\n";
}
/*
*------------------------------------------------------------------------
* Verify the norm calculation
*/
static void test_norms(const int vsize)
{
cout << "\n---> Verify norm calculations\n";
const double pattern = 10.25;
if( vsize % 2 == 1 )
_error("Sorry, size of the vector to test must be even for this test\n");
Vector v(vsize);
Vector v1(v);
cout << "\nAssign " << pattern << " to all the elements and check norms\n";
v = pattern;
cout << " 1. norm should be pattern*no_elems\n";
assert( v.norm_1() == pattern*v.q_no_elems() );
cout << " Square of the 2. norm has got to be pattern^2 * no_elems\n";
assert( v.norm_2_sqr() == sqr(pattern)*v.q_no_elems() );
cout << " Inf norm should be pattern itself\n";
assert( v.norm_inf() == pattern );
cout << " Scalar product of vector by itself is the sqr(2. vector norm)\n";
assert( v.norm_2_sqr() == v * v );
double ap_step = 1;
double ap_a0 = -pattern;
int n = v.q_no_elems();
cout << "\nAssign the arithm progression with 1. term " << ap_a0 <<
"\nand the difference " << ap_step << "\n";
register int i;
for(i=v.q_lwb(); i<=v.q_upb(); i++)
v(i) = (i-v.q_lwb())*ap_step + ap_a0;
int l = min(max((int)ceil(-ap_a0/ap_step),0),n);
double norm = (2*ap_a0 + (l+n-1)*ap_step)/2*(n-l) +
(-2*ap_a0-(l-1)*ap_step)/2*l;
cout << " 1. norm should be " << norm << "\n";
assert( v.norm_1() == norm );
norm = n*( sqr(ap_a0) + ap_a0*ap_step*(n-1) + sqr(ap_step)*(n-1)*(2*n-1)/6);
cout << " Square of the 2. norm has got to be "
"n*[ a0^2 + a0*q*(n-1) + q^2/6*(n-1)*(2n-1) ], or " << norm << "\n";
assert( v.norm_2_sqr() == norm );
norm = max(abs(v(v.q_lwb())),abs(v(v.q_upb())));
cout << " Inf norm should be max(abs(a0),abs(a0+(n-1)*q)), ie " << norm <<
"\n";
assert( v.norm_inf() == norm );
cout << " Scalar product of vector by itself is the sqr(2. vector norm)\n";
assert( v.norm_2_sqr() == v * v );
v1.clear();
compare(v,v1,"Compare the vector v with a zero vector");
cout << "\nConstruct v1 to be orthogonal to v as v(n), -v(n-1), v(n-2)...\n";
for(i=0; i<v1.q_no_elems(); i++)
v1(i+v1.q_lwb()) = v(v.q_upb()-i) * ( i % 2 == 1 ? -1 : 1 );
cout << "||v1|| has got to be equal ||v|| regardless of the norm def\n";
assert( v1.norm_1() == v.norm_1() );
assert( v1.norm_2_sqr() == v.norm_2_sqr() );
assert( v1.norm_inf() == v.norm_inf() );
cout << "But the scalar product has to be zero\n";
assert( v1 * v == 0 );
cout << "\nDone\n";
}
/*
*------------------------------------------------------------------------
* Root module
*/
main()
{
cout << "\n\n" << _Minuses <<
"\n\t\tVerify Operations on Vectors\n";
test_allocation();
test_element_op(20);
test_binary_op(20);
test_norms(20);
}